%**************************************************************************
%
%           "Endogenous Liquidity and Capital Reallocation"
%      Cui, W., Wright, R., & Zhu, Y. (2024), Journal of Political Economy
%                       Last Modified: Feb 2024.  
%
%**************************************************************************

%%%% Trend and business cycle properties of captial reallocation 

clear; close all; clc; format short

pics = 1; % show pictures

%% Step 1 ----------------------- Read data -------------------------------

time = 1971:1:2018;

%%%% NBER Recession data
Recessions = [ 
    1973+10/12,1975+2/12;
    1980,1980+6/12;
    1981+6/12,1982+10/12;
    1990+6/12,1991+2/12;
    2001+2/12,2001+10/12;
    2007+11/12,2009+5/12;];

%%%% 1971-2018 recession year index
Recession_ind = [0;0;1;1;1;0;0;0;0;1;1;1;0;0;0;0;0;0;0;1;1;...
    0;0;0;0;0;0;0;0;1;1;0;0;0;0;0;1;1;1;0;0;0;0;0;0;0;0;0];

%%%% Reallocation and Macro data
part_ratio          = importdata('participation.csv');              % from COMPUSTAT; our own calculations
deposit_rates       = importdata('deposit_rates_reallocation.csv'); % from call report; our own calculations
TFP_growthdata      = importdata('Fernald_TFP_growth.csv');         % from John Fernald website
reallocation_data   = importdata('reallocation_measures_Wei.csv');  % from COMPUSTAT;

load('FredData.mat') %%%% saved from before; you can run getDataOnline to see updated data

ID_87_18 = deposit_rates.data(:,1)<=2018;
transdep_rate_87_18 = deposit_rates.data(ID_87_18,2);

%%%% TFP
TFP_growth      = TFP_growthdata.data(:, 2) / 100;
A               = ones(length(TFP_growth),length(TFP_growth));
B               = tril(A);
log_TFP         = B * TFP_growth; 
log_TFP         = log_TFP - log_TFP(1);


%% Figure 1; Check R2 in the command Window
figure()
subplot(2,3,1)
hold on
x = log(inflation(18:end));
y = log((transdep_rate_87_18(1:end-1) + 1) ./ inflation(18:end));
scatter(x,y, 15, 'filled')
bestfit = polyfit(x,y,1);
yfit = polyval(bestfit, x);
plot(x,yfit)
ylabel('Transaction deposits')
theString = sprintf('Slope: %.2f', bestfit(1));
text(0.01,0,theString, 'FontSize', 14);

subplot(2,3,2)
hold on
x = log(inflation(2:end));
y = log((T_bill_1y(1:end-1) / 100 + 1) ./ inflation(2:end));
scatter(x,y, 15, 'filled')
bestfit = polyfit(x,y,1);
yfit = polyval(bestfit, x);
plot(x,yfit)
ylabel('1 Year T-Bill')
theString = sprintf('Slope: %.2f', bestfit(1));
text(0.08,0,theString, 'FontSize', 14);

subplot(2,3,3)
hold on
x = log(inflation(2:end));
y = log((T_bill_10y(1:end-1) / 100 + 1) ./ inflation(2:end));
scatter(x,y, 15, 'filled')
bestfit = polyfit(x,y,1);
yfit = polyval(bestfit, x);
plot(x,yfit)
ylabel('10 Year T-Bill')
theString = sprintf('Slope: %.2f', bestfit(1));
text(0.1,0,theString, 'FontSize', 14);

subplot(2,3,4)
hold on
x = log(inflation(2:end));
y = log((mortgage_rate(1:end-1) / 100 + 1) ./ inflation(2:end));
scatter(x,y, 15, 'filled')
bestfit = polyfit(x,y,1);
yfit = polyval(bestfit, x);
plot(x,yfit)
ylabel('30 Year Mortgage')
xlabel('Gross inflatioin (in log)','FontSize', 12)
theString = sprintf('Slope: %.2f', bestfit(1));
text(0.1,0.02,theString, 'FontSize', 14);

subplot(2,3,5)
hold on
x = log(inflation(2:end));
y = log((AAA_rate(1:end-1) / 100 + 1) ./ inflation(2:end));
scatter(x,y, 15, 'filled')
bestfit = polyfit(x,y,1);
yfit = polyval(bestfit, x);
plot(x,yfit)
ylabel('AAA')
xlabel('Gross inflatioin (in log)', 'FontSize', 12)
theString = sprintf('Slope: %.2f', bestfit(1));
text(0.1,0,theString, 'FontSize', 14);

subplot(2,3,6)
hold on
x = log(inflation(2:end));
y = log((BAA_rate(1:end-1) / 100 + 1) ./ inflation(2:end));
scatter(x,y, 15, 'filled')
bestfit = polyfit(x,y,1);
yfit = polyval(bestfit, x);
plot(x,yfit)
ylabel('BAA')
xlabel('Gross inflatioin (in log)', 'FontSize', 12)
theString = sprintf('Slope: %.2f', bestfit(1));
text(0.1,0.02,theString, 'FontSize', 14);



%% ------------------ Step 2: Data prepartion -----------------------------

%%%% Reallocation, capital expenditures, and debt Data 
%%%% (check against the .csv file)
R_share         = reallocation_data.data(:,2); % Wei's updated reallocation measures
P_share         = reallocation_data.data(:,3); % Wei's updated reallocation measures
debt            = reallocation_data.data(:,4);

disp('**********  Average R share and P share since 1984')
[mean(R_share(14:end)), mean(P_share(14:end))]
disp(' ')

%%%% participation data

part_ratio_data = part_ratio.data;

%%%% Nominal output data
C = Cons;
I = Inv + (G_inv - G_inv_def);
G = G_exp + G_inv_def;
Y = C + I + G ;

disp('**********  Nominal C-Y ratio, I-Y ratio, and G-Y between 1984 and 2018')
[mean(C(14:end)./Y(14:end)), mean(I(14:end)./Y(14:end)), mean(G(14:end)./Y(14:end))]
disp(' ')

ItoOutput = I(:) ./ (Y(:) - G(:));
CtoOutput = C(:) ./ (Y(:) - G(:));


%%%% real output data
G_real_I = G_real_E_plus_I - G_real_E - G_inv_def ./ G_inv_def_p * 100; % subtracting national defense
c        = C_real;
inv      = I_real + G_real_I;
output   = c + inv; % private output

%%%% money and debt data
id1             = part_ratio_data(:,1)>=1971 &  part_ratio_data(:,1)<=2018;
part_ratio_data = part_ratio_data(id1,2);

m2gdp           = money ./gdp;
debt_output     = debt ./gdp / 1000; %debt is in millions while output is in billions
debt            = debt ./ deflator_bv * 100;


%% -------------- Step 3: analyzing and plotting figure--------------------

%%%% Trend and Business cycles
BP_end_year         = 9;
starting_point      = 47; % 1971 is 47 numbers before 2018

% T_bill_1y_cycle     = bpass(log(T_bill_1y(end - starting_point: end)), 2, BP_end_year); %band pass filter 7-year business cycles: 2 years - 9 years
% T_bill_3m_cycle     = bpass(log(T_bill_3m(end - starting_point: end)), 2, BP_end_year);
% inflation_cycle     = bpass(log(inflation(end - starting_point: end)), 2, BP_end_year);
% R_share_cycle       = bpass(log(R_share(end - starting_point: end)), 2, BP_end_year);
% P_share_cycle       = bpass(log(P_share(end - starting_point: end)), 2, BP_end_year);
% debt_cycle          = bpass(log(debt(end - starting_point: end)), 2, BP_end_year);
% output_cycle        = bpass(log(output(end - starting_point: end)), 2, BP_end_year);
% c_cycle             = bpass(log(c(end - starting_point: end)), 2, BP_end_year);
% inv_cycle           = bpass(log(inv(end - starting_point: end)), 2, BP_end_year);
% Hours_cycle         = bpass(log(Hours(end - starting_point: end)), 2, BP_end_year);
% TFP_cycle           = bpass(log_TFP(end - starting_point: end), 2, BP_end_year);
% partr_cycle         = bpass(log(part_ratio_data(end - starting_point: end)), 2, BP_end_year);
% 
% T_bill_1y_trend     = log(T_bill_1y(end - starting_point: end)) - T_bill_1y_cycle;
% T_bill_3m_trend     = log(T_bill_3m(end - starting_point: end)) - T_bill_3m_cycle;
% inflation_trend     = log(inflation(end - starting_point: end)) - inflation_cycle;
% R_share_trend       = log(R_share(end - starting_point: end)) - R_share_cycle;
% P_share_trend       = log(P_share(end - starting_point: end)) - P_share_cycle;
% debt_trend          = log(debt(end - starting_point: end)) - debt_cycle;
% output_trend        = log(output(end - starting_point: end)) - output_cycle;
% partr_trend         = log(part_ratio_data(end - starting_point: end)) - partr_cycle;

T_bill_1y_cycle     = bpass(log(T_bill_1y(end - starting_point: end)), 2, BP_end_year); %band pass filter 6 quarters - 32 quarters, which is 2 years - 8 years
T_bill_3m_cycle     = bpass(log(T_bill_3m(end - starting_point: end)), 2, BP_end_year);
inflation_cycle     = bpass(log(inflation(end - starting_point: end)), 2, BP_end_year);
R_share_cycle       = bpass(log(R_share(end - starting_point: end)), 2, BP_end_year);
P_share_cycle       = bpass(log(P_share(end - starting_point: end)), 2, BP_end_year);
debt_cycle          = bpass(log(debt(end - starting_point: end)), 2, BP_end_year);
output_cycle        = bpass(log(output(end - starting_point: end)), 2, BP_end_year);
c_cycle             = bpass(log(c(end - starting_point: end)), 2, BP_end_year);
inv_cycle           = bpass(log(inv(end - starting_point: end)), 2, BP_end_year);
Hours_cycle         = bpass(log(Hours(end - starting_point: end)), 2, BP_end_year);
TFP_cycle           = bpass(log_TFP(end - starting_point: end), 2, BP_end_year);

ItoOutput_cycle     = bpass(ItoOutput(end - starting_point: end), 2, BP_end_year);
CtoOutput_cycle     = bpass(CtoOutput(end - starting_point:end), 2,  BP_end_year);

partr_cycle         = bpass(log(part_ratio_data(end - starting_point: end)), 2, BP_end_year);
debt_output_cycle   = bpass(log(debt_output(end - starting_point: end)), 2, BP_end_year);

T_bill_1y_trend     = log(T_bill_1y(end - starting_point: end)) - T_bill_1y_cycle;
T_bill_3m_trend     = log(T_bill_3m(end - starting_point: end)) - T_bill_3m_cycle;
inflation_trend     = log(inflation(end - starting_point: end)) - inflation_cycle;
R_share_trend       = log(R_share(end - starting_point: end)) - R_share_cycle;
P_share_trend       = log(P_share(end - starting_point: end)) - P_share_cycle;
debt_trend          = log(debt(end - starting_point: end)) - debt_cycle;
output_trend        = log(output(end - starting_point: end)) - output_cycle;
partr_trend         = log(part_ratio_data(end - starting_point: end)) - partr_cycle;
debt_output_trend   = debt_output-debt_output_cycle;
ItoOutput_trend     = ItoOutput - ItoOutput_cycle;
CtoOutput_trend     = CtoOutput - CtoOutput_cycle;

data = [R_share_cycle P_share_cycle debt_cycle T_bill_1y_cycle T_bill_3m_cycle, inflation_cycle];
csvwrite(['data_cycle_BP.txt'],data)

%%%% Computing statistics
G_C         = c_cycle;
G_output    = output_cycle;
G_I         = inv_cycle;
G_Hours     = Hours_cycle;
G_TFP       = TFP_cycle;
G_R_share   = R_share_cycle;
G_P_share   = P_share_cycle;
G_debt      = debt_cycle;
G_iota      = T_bill_1y_cycle;
G_pi        = inflation_cycle;
DtoOutput = debt./output;

display('********** output standard deviation')

std(G_output)
disp(' ')

display('********** Relative std.dev. (debt, R share, P share, Investment, and Inflation)')
[std(G_debt) / std(G_output),...
std(G_R_share)   / std(G_output),...
std(G_P_share)  / std(G_output),...
std(G_I) / std(G_output),...
std(G_pi) / std(G_output),...
]
disp(' ')

display('********** Correlation (debt, R share, P share, Investment, Inflation, Output)')
display('NaN means obmitted because of symmetry')
[1, corr(G_debt, G_R_share),...
corr(G_debt, G_P_share),...
corr(G_debt, G_I),...
corr(G_debt, G_pi),...
corr(G_debt, G_output);...
NaN, 1, corr(G_R_share, G_P_share),...
corr(G_R_share, G_I),...
corr(G_R_share, G_pi),...
corr(G_R_share, G_output);...
NaN, NaN, 1, corr(G_P_share, G_I),...
corr(G_P_share, G_pi),...
corr(G_P_share, G_output);...
NaN, NaN, NaN, 1, corr(G_I, G_pi),...
corr(G_I, G_output);...
NaN, NaN, NaN, NaN, 1, corr(G_pi, G_output);...
NaN, NaN, NaN, NaN, NaN, 1]
%corr(G_iota, G_Y)
disp(' ')

display('Correltaion of participation and output')
corr(output_cycle, partr_cycle)
disp(' ')

display('Correltaion of TFP and hours')
corr(G_Hours, G_TFP)
disp(' ')

display('******** Relative std and correlation (order: C, I, hours, TFP,  R-E ratio, P share, inflation)')
[std(G_C)    / std(G_output),...
std(G_I)    / std(G_output),...
std(G_Hours)/ std(G_output),...
std(G_TFP)  / std(G_output),...
std(G_R_share)   / std(G_output),...
std(G_P_share)  / std(G_output),...
std(G_pi) / std(G_output)]

[corr(G_C, G_output),...
corr(G_I, G_output),...
corr(G_Hours, G_output),...
corr(G_TFP, G_output),...
corr(G_R_share, G_output),...
corr(G_P_share, G_output),...
corr(G_pi, G_output)]

% save('Tbill3m_trend','T_bill_3m_trend')
% save('inflation','inflation')
% save('R_share','R_share')
% save('P_share','P_share')
% save('DtoOutput','DtoOutput');
% save('m2gdp','m2gdp');
% disp(' ')

%% Figures
if pics == 1
    %%%% Scatter Plots
    disp('********** Regression results')
    plotting_raw_trend_cycle_inflation 

    %%%% Time series Plotting
 
    %%%%
    figure()
    subplot(2,2,1)
    [AX,H1,H2] = plotyy(time, debt_cycle / std(debt_cycle), time, output_cycle / std(output_cycle));
    set(H1,'LineStyle','-','Linewidth',2,'color','k')
    set(H2,'LineStyle','--','Linewidth',2,'color','b')
    legend('Firm debt', 'output','Location','South')
    %ylabel('Reallocation / TotalInv', 'Leverage Index', 'Location','South')
    set(AX,'XLim',[1970 2020],'Xtick', [1970:5:2020], 'FontSize',12)
    set(AX(1),'Ycolor','k','Ylim',[-3,3],'Ytick',[-3:1:3])
    set(AX(1),'box','off')
    set(AX(2),'Ycolor','b','Ylim',[-3,3],'YTick',[-3:1:3])
    for r = 1:size(Recessions,1)
            ph(r) = patch([Recessions(r,1) Recessions(r,1) Recessions(r,2) Recessions(r,2)], ...
                          [get(gca,'YLim') fliplr(get(gca,'YLim'))], [0 0 0 0], 'k');
            set(ph(r),'FaceAlpha',0.2,'EdgeColor','none');
    end

    subplot(2,2,2)
    [AX,H1,H2] = plotyy(time, inv_cycle / std(inv_cycle), time, output_cycle / std(output_cycle));
    set(H1,'LineStyle','-','Linewidth',2,'color','r')
    set(H2,'LineStyle','--','Linewidth',2,'color','b')
    legend('New Investment', 'output','Location','south')
    %ylabel('Reallocation / TotalInv', 'Leverage Index', 'Location','South')
    set(AX,'XLim',[1970 2020],'Xtick', [1970:5:2020], 'FontSize',12)
    set(AX(1),'Ycolor','r','Ylim',[-3,3],'Ytick',[-3:1:3])
    set(AX(1),'box','off')
    set(AX(2),'Ycolor','b','Ylim',[-3,3],'YTick',[-3:1:3])
    for r = 1:size(Recessions,1)
            ph(r) = patch([Recessions(r,1) Recessions(r,1) Recessions(r,2) Recessions(r,2)], ...
                          [get(gca,'YLim') fliplr(get(gca,'YLim'))], [0 0 0 0], 'k');
            set(ph(r),'FaceAlpha',0.2,'EdgeColor','none');
    end

    subplot(2,2,3)
    [AX,H1,H2] = plotyy(time, R_share_cycle / std(R_share_cycle), time, output_cycle / std(output_cycle));
    set(H1,'LineStyle','-','Linewidth',2,'color','r')
    set(H2,'LineStyle','--','Linewidth',2,'color','b')
    legend('R share', 'output','Location','South')
    %ylabel('Reallocation / TotalInv', 'Leverage Index', 'Location','South')
    set(AX,'XLim',[1970 2020],'Xtick', [1970:5:2020], 'FontSize',12)
    set(AX(1),'Ycolor','r','Ylim',[-3,3],'Ytick',[-3:1:3])
    set(AX(1),'box','off')
    set(AX(2),'Ycolor','b','Ylim',[-3,3],'YTick',[-3:1:3])
    for r = 1:size(Recessions,1)
            ph(r) = patch([Recessions(r,1) Recessions(r,1) Recessions(r,2) Recessions(r,2)], ...
                          [get(gca,'YLim') fliplr(get(gca,'YLim'))], [0 0 0 0], 'k');
            set(ph(r),'FaceAlpha',0.2,'EdgeColor','none');
    end

    subplot(2,2,4)
    [AX,H1,H2] = plotyy(time, P_share_cycle / std(P_share_cycle), time, output_cycle / std(output_cycle));
    set(H1,'LineStyle','-','Linewidth',2,'color','m')
    set(H2,'LineStyle','--','Linewidth',2,'color','b')
    legend('P share', 'output','Location','South')
    %ylabel('Reallocation / TotalInv', 'Leverage Index', 'Location','South')
    set(AX,'XLim',[1970 2020],'Xtick', [1970:5:2020], 'FontSize',12)
    set(AX(1),'Ycolor','m','Ylim',[-3,3],'Ytick',[-3:1:3])
    set(AX(1),'box','off')
    set(AX(2),'Ycolor','b','Ylim',[-3,3],'YTick',[-3:1:3])
    for r = 1:size(Recessions,1)
            ph(r) = patch([Recessions(r,1) Recessions(r,1) Recessions(r,2) Recessions(r,2)], ...
                          [get(gca,'YLim') fliplr(get(gca,'YLim'))], [0 0 0 0], 'k');
            set(ph(r),'FaceAlpha',0.2,'EdgeColor','none');
    end

end



%% Compute medium-run and long-run

time = 1971:1:2018;

inflation_grid      = [0, 0.1/3, 0.2/3, 0.1]';

variables_trend     = [inflation_trend, ItoOutput_trend, CtoOutput_trend, R_share_trend, P_share_trend];
variables_itself    = [inflation, ItoOutput, CtoOutput, R_share, P_share];

save variables %%%% this saves everything
